home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / rays.i < prev    next >
Text File  |  1995-07-26  |  16KB  |  501 lines
  1. /*
  2.    HYPERB.I
  3.    Yorick functions to manipulate 3-D rays on a cylindrical mesh.
  4.  
  5.    $Id: rays.i,v 1.1 1993/09/02 22:48:38 munro Exp $
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /*
  11.    form_rays      -- convert from nrays lists to (5, 3, or 6)-by-nrays
  12.    best_rays      -- convert to 5-by-nrays representation
  13.    dirt_rays      -- convert to 3-by-nrays representation
  14.    internal_rays  -- convert to 6-by-nrays representation
  15.    picture_rays   -- return rays at centers of arbitrary "pixel mesh"
  16.  
  17.    get_s0         -- compute origins for slimits (internally, slimits
  18.                      is based on an s coordinate which is zero at the
  19.              point of closest approach to the origin, rather
  20.              than at the x,y,z point selected in the best_rays
  21.              or internal_rays coordinate systems)
  22.  
  23.    plray          -- plot a ray
  24.    plray_lims     -- set limits for plray
  25.  */
  26.  
  27. func form_rays(rays)
  28. /* DOCUMENT best= form_rays( [x, y, z, theta, phi] )
  29.          or dirt= form_rays( [x, y, theta] )
  30.      or internal= form_rays( [cos, sin, y, z, x, r] )
  31.      forms 5-by-nrays, 3-by-nrays, or 6-by-nrays ray representation
  32.      given individual lists of array coordinates.  The [...]
  33.      operator builds an nrays-by-5, nrays-by-3, or nrays-by-6
  34.      array, which form_rays transposes.  The "nrays" may represent
  35.      zero or more actual dimensions.
  36.    SEE ALSO: best_rays, dirt_rays, internal_rays, picture_rays
  37.  */
  38. { return transpose(rays, 2); }
  39.  
  40. func best_rays(rays)
  41. /* DOCUMENT best_rays(rays)
  42.      returns 5-element (x,y,z,theta,phi) representation of RAYS.
  43.      The first dimension of RAYS may be length 3, 5, or 6 to represent
  44.      the ray(s) in TDG/DIRT coordinates (x,y,theta), "best" coordinates
  45.      (x,y,z,theta,phi), or internal coordinates (cos,sin,y,z,x,r),
  46.      respectively.  The first dimension of the result always has length 5.
  47.  
  48.      The "best" coordinate system is the easiest to visualize:
  49.      (x,y,z) represents any point on the ray, while (theta,phi)
  50.      represents the ray direction in standard spherical coordinates
  51.      relative to the +z-axis.  Namely, theta is the angle from the
  52.      +z-direction to the ray direction (between 0 and pi), and phi is
  53.      the counterclockwise angle from the +x-axis to the projection of
  54.      the ray direction into the xy-plane, assuming xyz is a right-handed
  55.      coordinate system.
  56.  
  57.      As a specification of a ray, this system is doubly redundant because
  58.      the point (x,y,z) could be any point on the ray, and the underlying
  59.      mesh through which the ray propagates is cylindrically symmetric about
  60.      the z-axis.
  61.  
  62.      However, the slimits parameter -- used to specify the points along
  63.      a ray where the transport integration starts and stops -- is
  64.      measured from the point (x,y,z) specified as a part of the
  65.      (x,y,z,theta,phi) ray coordinate.  Thus, any change in the point
  66.      (x,y,z) on a ray must be accompanied by a corresponding change in
  67.      the slimits for that ray.
  68.  
  69.    SEE ALSO: form_rays, dirt_rays, internal_rays, get_s0, picture_rays
  70.  */
  71. {
  72.   dummy= use_origins(0);
  73.  
  74.   dim= dimsof(rays)(2);
  75.   if (dim==5) return rays;
  76.  
  77.   if (dim==3) {
  78.     /* convert from DIRT/TDG (x,y,theta) to (x,y,z,theta,phi) */
  79.     x= rays(1,..);
  80.     y= rays(2,..);
  81.     theta= rays(3,..);
  82.     return form_rays([x*cos(theta), y, x*sin(theta),
  83.               abs(theta), (theta>=0.0)*pi]);
  84.  
  85.   } else if (dim==6) {
  86.     /* convert from internal (cos,sin,y,z,x,r) to (x,y,z,theta,phi) */
  87.     theta= atan(rays(2,..), rays(1,..));
  88.     y= rays(3,..);
  89.     z= rays(4,..);
  90.     x= rays(5,..);
  91.     return form_rays([x, y, z, abs(theta), (theta<0)*pi]);
  92.   }
  93.  
  94.   return [];
  95. }
  96.  
  97. func dirt_rays(rays)
  98. /* DOCUMENT dirt_rays(rays)
  99.      returns 3-element (x,y,theta) representation of RAYS.
  100.      The first dimension of RAYS may be length 3, 5, or 6 to represent
  101.      the ray(s) in TDG/DIRT coordinates (x,y,theta), "best" coordinates
  102.      (x,y,z,theta,phi), or internal coordinates (cos,sin,y,z,x,r),
  103.      respectively.  The first dimension of the result always has length 3.
  104.  
  105.      The TDG/DIRT coordinate system is based on the coordinates (x,y)
  106.      in a plane normal to the ray.  Unfortunately, the old TDG and DIRT
  107.      codes used an angle theta which has the opposite sense from the
  108.      "best" and internal coordinates.  Therefore, conversion from
  109.      TDG/DIRT coordinates to internal coordinates will reverse the
  110.      sign of theta.  Conversion from TDG/DIRT coordinates to "best"
  111.      coordinates always results in positive theta, but the angle phi
  112.      will be pi for positive input theta.
  113.  
  114.      The slimits parameter -- used to specify the points along
  115.      a ray where the transport integration starts and stops -- is
  116.      measured from the point of closest approach of the ray described
  117.      by (x,y,theta) to the origin x=y=z=0.  Therefore, slimits is
  118.      independent of the TDG/DIRT ray coordinate representation.
  119.  
  120.    SEE ALSO: form_rays, best_rays, internal_rays, get_s0, picture_rays
  121.  */
  122. {
  123.   dummy= use_origins(0);
  124.  
  125.   dim= dimsof(rays)(2);
  126.   if (dim==3) return rays;
  127.  
  128.   if (dim==5) {
  129.     /* convert from best (x,y,z,theta,phi) to DIRT/TDG (x,y,theta) */
  130.     x= rays(1,..);
  131.     y= rays(2,..);
  132.     z= rays(3,..);
  133.     theta= rays(4,..);
  134.     phi= rays(5,..);
  135.     cosp= cos(phi);
  136.     sinp= sin(phi);
  137.     sgn= double((cosp>=0.0)*2-1);
  138.     return form_rays(sgn*[(x*cosp+y*sinp)*cos(theta) - z*sin(theta),
  139.               (y*cosp-x*sinp), -theta]);
  140.  
  141.   } else if (dim==6) {
  142.     /* convert from internal (cos,sin,y,z,x,r) to DIRT/TDG (x,y,theta) */
  143.     cost= rays(1,..);
  144.     sint= rays(2,..)
  145.     y= rays(3,..);
  146.     z= rays(4,..);
  147.     x= rays(5,..);
  148.     /* NOTE SIGN CHANGE IN THETA */
  149.     return form_rays([x*cost-z*sint, y, -atan(sint, cost)]);  
  150.   }
  151.  
  152.   return [];
  153. }
  154.  
  155. func internal_rays(rays)
  156. /* DOCUMENT internal_rays(rays)
  157.      returns 6-element (cos,sin,y,z,x,r) representation of RAYS.
  158.      The first dimension of RAYS may be length 3, 5, or 6 to represent
  159.      the ray(s) in TDG/DIRT coordinates (x,y,theta), "best" coordinates
  160.      (x,y,z,theta,phi), or internal coordinates (cos,sin,y,z,x,r),
  161.      respectively.  The first dimension of the result always has length 6.
  162.  
  163.      The internal coordinates are what Drat uses internally to
  164.      describe the ray.  The coordinate system is rotated about the
  165.      z-axis until the ray lies in a plane of constant y (there are at
  166.      least two ways to do this).  The point (x,y,z) can be any point on
  167.      the ray, and r=sqrt(x^2+y^2) is the corresponding cylindrical radius.
  168.      The clockwise angle theta from the +z-axis to the ray direction
  169.      (which always lies in the zx-plane) determines cos=cos(theta) and
  170.      sin=sin(theta).
  171.  
  172.      As a specification of a ray, this system is triply redundant because
  173.      the point (x,y,z) could be any point on the ray, both the sine and
  174.      cosine of theta appear, and r=sqrt(x^2+y^2).
  175.  
  176.      However, the slimits parameter -- used to specify the points along
  177.      a ray where the transport integration starts and stops -- is
  178.      measured from the point (x,y,z) specified as a part of the
  179.      (cos,sin,y,z,x,r) ray coordinate.  Thus, any change in the point
  180.      (x,y,z) on a ray must be accompanied by a corresponding change in
  181.      the slimits for that ray.
  182.  
  183.    SEE ALSO: form_rays, best_rays, dirt_rays, get_s0, picture_rays
  184.  */
  185. {
  186.   dummy= use_origins(0);
  187.  
  188.   dim= dimsof(rays)(2);
  189.  
  190.   if (dim==6) {
  191.     /* assure that r is consistent with x and y */
  192.     rays(6,..)= abs(rays(3,..), rays(5,..));
  193.     return rays;
  194.  
  195.   } else if (dim==5) {
  196.     /* convert from best (x,y,z,theta,phi) to internal (cos,sin,y,z,x,r) */
  197.     x= rays(1,..);
  198.     y= rays(2,..);
  199.     z= rays(3,..);
  200.     theta= rays(4,..);
  201.     phi= rays(5,..);
  202.     cosp= cos(phi);
  203.     sinp= sin(phi);
  204.     sgn= double((cosp>=0.0)*2-1);
  205.     return form_rays([cos(theta), sgn*sin(theta), sgn*(y*cosp-x*sinp), z,
  206.               sgn*(x*cosp+y*sinp), abs(x, y)]);
  207.  
  208.   } else if (dim==3) {
  209.     /* convert from DIRT/TDG (x,y,theta) to internal (cos,sin,y,z,x,r) */
  210.     x= rays(1,..);
  211.     y= rays(2,..);
  212.     theta= rays(3,..);
  213.     cost= cos(theta);
  214.     sint= sin(theta);
  215.     xcost= x*cost;
  216.     /* NOTE SIGN CHANGE IN THETA */
  217.     return form_rays([cost, -sint, y, x*sint, xcost, abs(xcost, y)]);
  218.   }
  219.  
  220.   return [];
  221. }
  222.  
  223. func get_s0(rays)
  224. /* DOCUMENT get_s0(rays)
  225.      returns the s-coordinate of the point of closest approach of
  226.      the RAYS to the origin x=y=z=0.  The length of the first dimension
  227.      of RAYS may be either 3, 5, or 6; this first dimension will not
  228.      be present in the result.
  229.  
  230.      The s-coordinate represents distance along the ray, increasing in
  231.      the direction the ray moves.  The 5 and 6 component ray coordinates
  232.      include a reference point (x,y,z) on the ray; s=0 at that point.
  233.      For the 3 component ray coordinate, get_s0 always returns 0.
  234.  
  235.    SEE ALSO: best_rays, dirt_rays, internal_rays
  236.  */
  237. {
  238.   dummy= use_origins(0);
  239.  
  240.   dim= dimsof(rays)(2);
  241.  
  242.   if (dim==3) {
  243.     return array(0.0, dimsof(rays(1,..)));
  244.  
  245.   } else if (dim==5) {
  246.     x= rays(1,..);
  247.     y= rays(2,..);
  248.     z= rays(3,..);
  249.     theta= rays(4,..);
  250.     phi= rays(5,..);
  251.     cost= cos(theta);
  252.     sint= sin(theta);
  253.     x= x*cos(phi)+y*sin(phi);
  254.  
  255.   } else if (dim==6) {
  256.     cost= rays(1,..);
  257.     sint= rays(2,..)
  258.     z= rays(4,..);
  259.     x= rays(5,..);
  260.   }
  261.  
  262.   return -z*cost-x*sint;
  263. }
  264.  
  265. func picture_rays(xpict, ypict, ray, theta_up, phi_up)
  266. /* DOCUMENT picture_rays(xpict, ypict, ray)
  267.          or picture_rays(xpict, ypict, ray, theta_up, phi_up)
  268.      returns 2-D array of rays, one at the center of each zone (which
  269.      represents a pixel here) of the mesh (XPICT, YPICT).  The rays are
  270.      all parallel to the given RAY (a 3, 5, or 6 element vector).  The
  271.      (XPICT, YPICT) coordinates are in the plane perpendicular to the rays,
  272.      with the origin XPICT=YPICT=0 at the given RAY.
  273.  
  274.      If (THETA_UP, PHI_UP) are given, the +YPICT-axis will lie along the
  275.      projection of the (THETA_UP, PHI_UP) direction into the (XPICT, YPICT)
  276.      plane.  The default (THETA_UP, PHI_UP) is (pi/2, pi/2) -- the +y-axis
  277.      -- unless (THETA, PHI) is the y-axis, in which case it is (pi/2, 0)
  278.      -- the +x-axis.  This matches the DIRT/TDG ray coordinate convention
  279.      in the sense that if RAY is [0,0,theta], then
  280.      (zncen(XPICT),zncen(YPICT)) are the DIRT/TDG (x,y) coordinates for
  281.      the rays.
  282.  
  283.      If XPICT and YPICT are imax-by-jmax, the returned array will have
  284.      dimensions 5-by-(imax-1)-by-(jmax-1).  That is, "best" coordinates
  285.      are returned.  The (x,y,z) of all of the returned rays will lie in
  286.      the plane perpendicular to the ray passing through the given central
  287.      RAY.
  288.  
  289.    SEE ALSO: form_rays, best_rays, dirt_rays, internal_rays, area
  290.  */
  291. {
  292.   dummy= use_origins(0);
  293.   ray= best_rays(ray);
  294.   x= ray(1);
  295.   y= ray(2);
  296.   z= ray(3);
  297.   theta= ray(4);
  298.   phi= ray(5);
  299.  
  300.   /* form (rx,ry,rz) ray direction cosines */
  301.   sint= sin(theta);
  302.   rx= sint*cos(phi);
  303.   ry= sint*sin(phi);
  304.   rz= cos(theta);
  305.  
  306.   /* form (ux,uy,uz) picture up direction cosines */
  307.   if (!is_void(theta_up)) {
  308.     sint= sin(theta_up);
  309.     ux= sint*cos(phi_up);
  310.     uy= sint*sin(phi_up);
  311.     uz= cos(theta_up);
  312.   } else {
  313.     if (rz || rx) { ux= 0.0; uy= 1.0; }
  314.     else { ux= 1.0; uy= 0.0; }
  315.     uz= 0.0;
  316.   }
  317.  
  318.   /* form (px,py,pz) picture right direction cosines */
  319.   px= uy*rz - uz*ry;
  320.   py= uz*rx - ux*rz;
  321.   pz= ux*ry - uy*rx;
  322.   len= abs(px, py, pz);
  323.   px/= len;
  324.   py/= len;
  325.   pz/= len;
  326.  
  327.   /* project (ux,uy,uz) perpendicular to (rx,ry,rz) */
  328.   ux= ry*pz - rz*py;
  329.   uy= rz*px - rx*pz;
  330.   uz= rx*py - ry*px;
  331.  
  332.   /* get coordinates of zone centers in picture plane */
  333.   xp= xpict(zcen, zcen);
  334.   yp= ypict(zcen, zcen);
  335.  
  336.   return form_rays([x+xp*px+yp*ux, y+xp*py+yp*uy, z+xp*pz+yp*uz,
  337.             theta, phi]);
  338. }
  339.  
  340. func plray(ray)
  341. /* DOCUMENT plray, ray
  342.      where RAY is a vector of length 5 representing [x, y, z, theta, phi]
  343.      -- a point (x,y,z) on the ray, and the ray direction (theta,phi)
  344.      relative to the z-axis.  The ray hyperbola is plotted with the
  345.      z-axis horizontal.  The portion of the ray plotted is determined
  346.      by the plray_lims command, which must be issued prior to the first
  347.      plray.
  348.      The 3 and 6 component ray formats are also accepted.
  349.    SEE ALSO: plray_lims, drat_rays, dirt_rays, internal_rays
  350.  */
  351. {
  352.   dims= dimsof(ray);
  353.   if (dims(1)!=1) error, "ray must be a vector of length 5 (or 3 or 6)";
  354.   if (dims(2)!=5) ry= best_rays(ray);
  355.   else ry= ray;
  356.  
  357.   dummy= use_origins(0);
  358.  
  359.   x0= ry(1);
  360.   y0= ry(2);
  361.   z0= ry(3);
  362.   theta= ry(4);
  363.   phi= ry(5);
  364.  
  365.   sint= sin(theta);
  366.   cost= cos(theta);
  367.   sinp= sin(phi);
  368.   cosp= cos(phi);
  369.  
  370.   rtp= abs(x0*sinp - y0*cosp);  /* turning point radius */
  371.  
  372.   /**/ extern plray_zx, plray_zn, plray_rx;
  373.   /**/ extern plray_n;   /* always odd! */
  374.  
  375.   if (sint < 1.e-4) {
  376.     /* rays which are nearly parallel to the axis need special treatment */
  377.     z= span(plray_zn, plray_zx, 5);
  378.     if (cost<0.0) z= z(::-1);
  379.     r= (z-z0)*sint/cost;
  380.     r= abs(x0+r*cosp, y0+r*sinp);
  381.  
  382.   } else {
  383.     /* form th == tan(angle)/abs(tan(theta)) where angle is evenly
  384.        spaced as th varies from -1 to 1 */
  385.     abstan= abs(sint)/(abs(cost)+1.e-20);
  386.     absth= atan(abstan);
  387.     fuzz= 1.e-6;
  388.     th= span(-absth, absth, plray_n);
  389.     th(1)= -(1.0-0.5*fuzz*fuzz);
  390.     th(2:-1)= tan(th(2:-1))/abstan;
  391.     th(0)= -th(1);
  392.  
  393.     sqth= sqrt((1.0-th)*(1.0+th));
  394.     sqth(1)= sqth(0)= fuzz;
  395.  
  396.     /* rays which nearly hit the axis need special treatment */
  397.     if (rtp < fuzz*plray_rx) rt= fuzz*plray_rx/sqth;
  398.     else rt= rtp/sqth;
  399.  
  400.     zt= z0 + (rt*th - (x0*cosp-y0*sinp)*sign(sint))*sign(cost)/abstan;
  401.  
  402.     /* clip (zt,rt) to (z,r) with z between plray_zn and plray_zx */
  403.     list= where(zt>plray_zn & zt<plray_zx);
  404.     if (min(abs(zt(dif))) > 0.0) {
  405.       if (cost>=0.0) {  /* z increasing */
  406.         if (min(zt)<=plray_zn) {
  407.       z= [plray_zn];
  408.       r= interp(rt, zt, z);
  409.     } else {
  410.       z= [];
  411.       r= [];
  412.     }
  413.     grow, z, zt(list);
  414.     grow, r, rt(list);
  415.         if (max(zt)>=plray_zx) {
  416.       grow, z, [plray_zx];
  417.       grow, r, interp(rt, zt, [plray_zx]);
  418.     }
  419.       } else {          /* z decreasing */
  420.         if (max(zt)>=plray_zx) {
  421.       z= [plray_zx];
  422.       r= interp(rt, zt, z);
  423.     } else {
  424.       z= [];
  425.       r= [];
  426.     }
  427.     grow, z, zt(list);
  428.     grow, r, rt(list);
  429.         if (min(zt)<=plray_zn) {
  430.       grow, z, [plray_zn];
  431.       grow, r, interp(rt, zt, [plray_zn]);
  432.     }
  433.       }
  434.     } else {
  435.       z= zt(list);
  436.       r= rt(list);
  437.     }
  438.   }
  439.  
  440.   /* clip (z,r) to r < rx */
  441.   if (numberof(r) && numberof((list= where(r<plray_rx)))) {
  442.     if (numberof(list) < numberof(r)) {
  443.       rt= r;
  444.       zt= z;
  445.       i= rt(mnx);  /* the interp operation seems to be safe even if many
  446.               of the rt are clustered at this minimum value */
  447.       if (min(list)<=i && max(rt(:i))>=plray_rx) {
  448.     r= [plray_rx];
  449.     z= interp(zt(:i), rt(:i), r);
  450.       } else {
  451.     r= [];
  452.     z= [];
  453.       }
  454.       grow, r, rt(list);
  455.       grow, z, zt(list);
  456.       if (max(list)>=i && max(rt(i:))>=plray_rx) {
  457.     grow, r, [plray_rx];
  458.     grow, z, interp(zt(i:), rt(i:), [plray_rx]);
  459.       }
  460.     }
  461.  
  462.     /* finally, plot the result */
  463.     /* Note: Gist substitutes the marker for \001 if that is the
  464.              first character of the legend.  */
  465.     legend=
  466.       swrite(format="\001: theta= %.4g, phi= %.4g, thru (%.4g, %.4g, %.4g)",
  467.          theta, phi, x0, y0, z0)(1);
  468.     plg, r, z, rays=1, legend=legend;
  469.  
  470.   } else {
  471.     write, "WARNING (plray) ray not within limits (use plray_lims)";
  472.   }
  473. }
  474.  
  475. func plray_lims(zmin, zmax, rmax, rx)
  476. /* DOCUMENT plray_lims, zmin, zmax, rmax
  477.      sets the (z,r) coordinate limits for the plray command.  Subsequent
  478.      rays will be clipped to the region from z=ZMIN to z=ZMAX, and from
  479.      r=0 to r=RMAX.
  480.    SEE ALSO: plray
  481.  */
  482. {
  483.   /**/ extern plray_zx, plray_zn, plray_rx;
  484.   if (!is_void(rx)) rmax= rx;  /* allow for limits syntax as well */
  485.   if (zmin!=zmax && rmax>0.0) {
  486.     plray_zx= double(max(zmin, zmax));
  487.     plray_zn= double(min(zmin, zmax));
  488.     plray_rx= double(rmax);
  489.     limits, plray_zn, plray_zx, 0.0, plray_rx;
  490.   } else {
  491.     error, "zmin==zmax or rmax<=0.0 are not legal ray limits";
  492.   }
  493. }
  494.  
  495. plray_zx= 1.0;
  496. plray_zn= -1.0;
  497. plray_rx= 1.0;
  498. plray_n= 31;    /* This should always be an odd number to ensure that
  499.            the turning point of the ray is always one of the
  500.            points plotted.  */
  501.